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The electronic structure near defects (such as impurities) in superconductors is explored using a 
new, fully self-consistent technique. This technique exploits the short-range nature of the impurity 
potential and the induced change in the superconducting order parameter to calculate features in 
the electronic structure down to the atomic scale with unprecedented spectral resolution. Magnetic 
and non-magnetic static impurity potentials are considered, as well as local alterations in the pairing 
interaction. Extensions to strong-coupling superconductors and superconductors with anisotropic 
order parameters are formulated. 



I. INTRODUCTION 

Low-tcniperature superconductors almost always have 
a high concentration of non-magnetic impurities. Even 
in the dirty limit, however, where the mean free path 
is shorter than the coherence length, superconductivity 
endures"'^ . This phenomenon can be understood by gener- 
alizing BCS pairing^ to be between degenerate Kramers 
partnered states in a time-reversal invariant system'^. 
Magnetic impurities, which do break time- reversal invari- 
ance, have more profound effects on the superconductor 
in dilute concentrations than non-magnetic impurities, 
lowering the critical temperature Tc^^^ and producing 
localized states within the gap^~^^ which at sufficient 
concentrations hybridize to produce gaplessness^'*. In 
the course of the investigation of the effects of impuri- 
ties on superconductivity during the last four decades, 
the primary emphasis has been the influence of impuri- 
ties on bulk properties. These effects have been treated 
within a strong-coupling formalism (e.g. Ref. 15), but 
only very recently self-consistently and beyond the Born 
approximation^^. The above work was primarily con- 
cerned with bulk or impurity-averaged characteristics 
and ignored the spatial structure of electronic properties 
near to the impurity. 

Among the first local properties calculated in the vicin- 
ity of an impurity in a superconductor were the structures 
of screening clouds around a charged impurity^^'^^ and 
a magnetic impurity^*'^^ in a superconductor (charac- 
terized by exponentially-decaying Friedel-like^" oscilla- 
tions). The oscillation of the order parameter around 
a magnetic impurity was first evaluated without self- 
consistency^^"^'^. A self-consistent calculation of the or- 
der parameter at the impurity and very far away for weak 
impurity potentials was reported by Schlottmann^'^ for 
magnetic impurities and Fetter^^ for non-magnetic im- 
purities. 

Interest in local properties near impurities in super- 
conductors has been revived by the capability of scan- 
ning tunneling microscopy (STM) to perform localized 



spectroscopic measurements. The spatially-dependent 
differential conductivity around a non-magnetic impu- 
rity at the surface of a superconductor has been con- 
sidered theoretically^^ for both isotropic and anisotropic 
order parameters. This differential conductivity is pro- 
portional to the local density of states (LDOS) around 
the impurity. In Ref. 26 the impurity was modeled 
as a point defect, and spatial oscillations in the LDOS 
at various voltages were calculated. These oscillations 
are the superconductor's analog of oscillations in the 
LDOS created by step edges and defects on noble metal 
surfaces^^'^*. The superconductor's LDOS oscillations 
would allow one to measure the anisotropy of the su- 
perconductor's order parameter. The conductance oscil- 
lations for voltages just above a gap minimum or maxi- 
mum are strongly pronounced in the real-space directions 
corresponding to the momenta of the gap minimum and 
maximum. Calculations followed which considered sharp 
energy features in the scattering process, such as reso- 
nant states^^'^°. Nevertheless, an important assumption 
of these calculations has remained unchecked in detail, 
that the electronic distortions induced by the impurity 
are local, including the deformation of the order parame- 
ter. Self-consistent calculations using the Bogohubov-de 
Gennes (BdG) equations^^ followed for two-dimensional 
systems'^^"'^^ , but have been hampered by finite-size ef- 
fects. 

A magnetic impurity differs from a non-magnetic im- 
purity in that a localized state exists around it^"^'^''^^. 
The first calculations of the LDOS of the localized and 
continuum states around a magnetic impurity were per- 
formed recently both with a simplified analytic model 
and numerically via the new self-consistent technique'^^ 
described in this paper. Calculations of the LDOS of the 
localized state with angular momentum quantum num- 
ber £ — were performed with a slightly different ana- 
lytic model and numerically in two dimensions via the 
self-consistent BdG equations'^^. These calculations were 
motivated by preliminary experimental results around a 
Gd atom on a niobium surface^^. 
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Another type of defect is a local change in the pair- 
ing interaction. The resulting spatially-dependent order 
parameter can then distort the LDOS. An order param- 
eter suppression can even localize states, as in a vortex 
core^O'^i. 

The local electronic properties of all of these de- 
fects can be calculated self-consistently from the Gor'kov 
equation^^ without further approximation with the new 
technique introduced in Ref. 37. Our technique for calcu- 
lating the electronic structure around a defect in a super- 
conductor is related to the Koster-Slater inversion tech- 
niques for determining the local electronic structure of 
impurities in metals**^'^^. Since its original application 
to localized vibrational modes*^, this algorithm has been 
applied to numerous problems including deep levels in 
semiconductors'*® and impurity states in magnets^''. The 
new Koster-Slater technique separates space around the 
defect into two regions: the near field and the far field. 
The far field is a region distant enough from the defect 
that the potential is insignificant and the order parame- 
ter has relaxed back to its homogeneous value. The near 
field is the region close to the defect where the potential is 
finite or the order parameter is distorted. In essence, the 
Gor'kov equation that determines the Green's functions 
of the inhomogeneous superconductor is inverted in the 
real space region of the near field. This paper describes 
the technique in detail, expands on an analytic model 
introduced in Ref. 37 and reports several calculations of 
the properties of the defects described above. 

Previous attempts to calculate the structure of the near 
field have often used the BdG equations'"*^. Other theories 
of inhomogeneities, such as Ginzburg-Landau theory'*^ or 
the Eilenberger equations^®, treat the spatial degrees of 
freedom as coarse-grained over the superconductor's co- 
herence length. Coarse-grained approximations are not 
appropriate for considering electronic structure on the 
atomic scale near a defect. The BdG equations are gen- 
eralized Schrodinger equations for the electron and hole 
wavefunctions of a quasiparticle, and are valid for a su- 
perconductor with an arbitrarily-varying order parame- 
ter, only constrained by the validity of BCS theory. 

Unfortunately, these equations have significant prac- 
tical difficulties as well. Despite qualitative success 
modeling STM measurements of a single vortex in su- 
perconducting NbSe2^°, calculations of the electronic 
structure^^ using the BdG equations are hampered by 
the difference in energy scales between the Fermi energy 
and the order parameter. Since the BdG equations are 
solved numerically for a finite system, the difficulty of 
the calculation is determined by the necessary spectral 
resolution. The key energy scale which must be resolved 
is the superconducting gap. Thus, the numerical diffi- 
culty increases as the ratio of the Fermi energy to the 
gap becomes large. Hence the band structure assumed 
for the superconductor must be somewhat unrealistic (for 
Rcfs. 54, 55 the Fermi wavelength was approximately 
lOOA, which is inappropriately large^®). This limitation 
extends to calculations of the interaction between a vor- 



tex and an impurity^^, the characteristics of the vortex 
lattice^*'^^, and work on a non-magnetic impurity^^"^^ 
and a magnetic impurity^^ in a two-dimensional s-wave 
or d-w&ve superconductor. In contrast, the computa- 
tional requirements of the Koster-Slater technique are 
determined by the range of the impurity potential, rather 
than the necessary spectral precision. 

In the Section II of this paper we first describe the 
BdG formalism for local defect potentials and then com- 
pare with the new, Koster-Slater formalism. Section III 
describes an analytic model, based on a delta-function 
potential, which reproduces some of the quantitative be- 
havior of the numerical results. Section IV discusses the 
results of the numerical calculations for magnetic impu- 
rities, non-magnetic impurities, impurities incorporating 
both magnetic and non-magnetic potentials, and inhomo- 
geneities in the pairing interaction. A heuristic picture of 
the electronic structure near these impurities will be pre- 
sented here, and the calculations will be compared with 
the analytic model of Section III. Section V generalizes 
the formalism of Section II to the case of strong-coupling 
and anisotropically-paired superconductors. 



II. FORMALISM 

A. Bogoliubov-de Gennes Equations 

To place our new formalism in context we will contrast 
it with the BdG equations, which are Schrodinger-like 
equations for the electron and hole components of the 
quasiparticle wavefunctions w(x) and v{x.) respectively. 
These are, for a free electron band structure with mass 
m, the positive-energy (£■) solutions to 



2m 
2m 



E + Vo{^)+aVs{^) 
E-Vo{^)+aVs{^) 



u<,(x) -h A(x)v^(x)= 

w^(x) A(x)u<,(x)= 0. 

(1) 



Here aVs is a position-dependent, spin-dependent po- 
tential, such as one originating from an impurity with a 
classical spin. Vq is a position-dependent non-magnetic 
potential and A(x) is the inhomogeneous order parame- 
ter. A(x) can be chosen real since the defect potential is 
real®°. The quantization direction of the electronic spins 
in the superconductor (ct = ±1/2) is chosen parallel to 
the classical spin. A classical spin has no quantum dy- 
namics, and cannot flip the quasiparticle spin. Hence 
spin is a good quantum number for the cjuasiparticles 
and only two coupled equations (Eqs. (1)) are required. 
The combinations aVs ± Vq have physical significance: 
o'Vs -|- Vo is the potential felt by an electron of spin a, 
while aVs — Vq is the potential felt by a hole of spin a. 

The spatially-dependent order parameter is deter- 
mined self-consistently: 
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A(x) = 5^ 7(x, i;„<.)u„<.(x)<,(x)tanh ( ^ j (2) 

no- \ B / 

where n labels the states for each spin a, T is the temper- 
ature, and fc_B is Boltzmann's constant. ^{x,Ena) is the 
effective electron-electron interaction potential, which is 

7(x, Ena)= lo Ena < TVjJD 

= Ena> huJD- (3) 

For a spherically-symmetric defect the wavefunctions 
are eigenstates of the angular momentum with quantum 
numbers £ and m. Typically the defect is placed in a 
sphere of radius R with appropriate boundary conditions. 
The value of R is determined by the spectral resolution 
necessary for accurately evaluating Eq. (2) and the spec- 
tral width of features measurable by (for example) the 
STM. The typical complications resulting from approx- 
imating an infinite system by a finite-size system apply, 
such as discrete states above the energy gap and the 
heavy investment of computer time required for large 
values of R. For example, in the calculations for the 
vortex in NbSe2^^'^^, cf/Aq = 32 was the largest ratio 
of the Fermi energy to the homogeneous order parame- 
ter considered. This value is unrealistic, and is a result 
of inappropriately fitting the coherence length and up- 
per critical field of NbSe2 to a free-electron model. A 
more realistic band structure^^ has a bandwidth to en- 
ergy gap ratio at least an order of magnitude greater, 
and a highly anisotropic Fermi velocity. A system where 
the frec-clcctron model is appropriate, such as niobium, 
has an ef/Aq = 705 and is numerically inaccessible for 
general potentials Vs and Vq. 

B. Koster-Slater Inversion Formalism 

We now introduce a self-consistent method which 

works within a sphere whose radius is determined by the 
range of the defect's potential and utilizes the contin- 
uum spectrum of the homogeneous superconductor. In 
essence, we invert the Gor'kov equation in real space. 
The Gor'kov equation^^ for a defect in a superconductor 
can be written in the Nambu formalism^^ as: 

j dx" [(5(x - x") - g(x, x"; c^)V(x")] G(x", x; u;) 

= g(x,x';w) (4) 
where the inhomogeneous retarded Green's function, 

G(x,x';a;)=f^l("'";'"j ^i"'"''") (5) 
^ ' ' ^ \ (x,x ,(jj) — GJ(x,x ,— w) J ^ ' 

The elements of this matrix are 

G|(x,x';a;) = 

/CXD 
die^-*^(i)(0|{VT (x'; t), V| (x; 0)}|0), (6) 
-CXD 



i^(x,x';a;) = 

/OO 
dte^-*^(i)(0|{VT(x'; i). V'i(x; 0)}|0), (7) 
-OO 

F*(x,x';a;) = 

/OO 
rfte'"*0(t)(O|{^|(x'; t), V| (x; 0)}|0), (8) 
-OO 

-GJ(x,x';-u;) = 

/OO 
dte^-*^(i)(0|{Vl (x'; i), V'i(x; 0)}|0). (9) 
-oo 

The explicit subscripts \ and j do not refer to the spin 
of the excitation in the superconductor but rather to the 
spin band of the normal state used to construct the exci- 
tation. The key concept is that the spin-up band contains 
both u'p electrons and down holes just as the spin- down 
band contains both down electrons and wp holes. The 
convention employed here is standard in the theory of 
semiconductors where a spin-up electron excited above 
the Fermi energy leaves a spin-down hole below the Fermi 
energy. This is convenient for magnetic potentials since 
if spin-up electrons are attracted to a magnetic impurity 
spin-down holes should be repelled by the impurity. In 
the presence of a single classical impurity spin, the quasi- 
particle spin is a good quantum number despite electron- 
hole mixing. Our convention determines the composition 
of a spin-up quasiparticlc to be part spin-up electron and 
part spin-up hole, rather than part spin-up electron and 
part spin-down hole. 

For a; > 0, G(x,x':i.j) describes spin-up excitations, 
involving the mixing of electrons in the spin-up band 
with holes in the spin-down band. For lo < 0, G(x, x'; w) 
describes spin-down excitations, involving the mixing of 
electrons in the spin-down band with holes in the spin-up 
band. Since spin is a good quantum number it is not nec- 
essary to use a 4 X 4 formalism, such as that of Ref. 62. 
The notation here has been simplified and improved rel- 
ative to Ref. 37. 

The homogeneous Green's function g is independent 
of a, so 

g(x, x» = f fif ' ^ ; ^\). (10) 

&\ , , J \^ /*(x,x';a;) -g*{x,x';-uj) J ^ ' 
The potential 

'^y^'l-y ^A(x") ys(x") - Vo(x") j ^^^> 

where (5A(x) = A(x) — Aq and Vs, Vq, and A(x) have 
similar meaning above as in the BdG equations. The 
factor of one-half from the electron spin has been incor- 
porated into the potential Vs- To be concrete, and with- 
out loss of generality, the spin direction attracted by the 
potential will be called spin up, and the spin direction 
repelled will be called spin down^^. The self-consistency 
equation for the order parameter is 
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/oo poo / 

c?a;7(x,a;)n(u;)ImF(x, x;a;) — (12) 5(x) = / dujn{uj) j - 
oo J — oo \ 

_ J -oo \ 



Im [G^ (x, x; w) - Gi (x, x; w)] ' 



where n{uj) is the Fermi occupation function at temper- 
ature T. 

Self-consistent spin-dependent and charge-dependent 
potentials can also be constructed for V(x) using the 
calculated spatial structure of the spin s(x) and charge 
p(x) around the defect, 



ImGo-(x, x;a;) 



(13) 
(14) 



No calculations are reported for such potentials'^ in this 
paper. 



J 



For a BCS superconductor with an isotropic gap in a parabolic band, for w much smaller than the Fermi energy, 

_e V"o '/T^oc. coskFr + 



/(x, x';uj) 



kpr 



: sin fcirr 



(15) 



where r = |x — x'| and No is the density of states for 
each spin at the Fermi level. The coherence length, 
^ = hvp /t^^o, where vp \s the Fermi velocity. These 
expressions are valid for uj above and below so long 
as the imaginary parts of both / and g are multiplied by 
sgnw. 

One strength of this formalism is its reliance on the 
short-range nature of the defect's potential. Solution of 
Eq. (4) requires inverting the frequency-dependent real- 
space matrix 

M(x, x'; w) = <5(x - x') - g(x, x'; a;)V(x'). (16) 

The structure of M(x, x'; w) allows for a precise descrip- 
tion of the difference between the near field and the far 
field. We require that the defect's potential V(x) is zero 
for |x| > R. The space |x| > R belongs to the far field, 
whereas the space |x| < R belongs to the near field. We 
can then separate any real-space matrix A symbolically 
into four pieces: 



J 



(17) 



where n and / label the near-field region and far-field re- 
gion, respectively. The particular example of M is block- 
triangular: 



M 



I - g"^"V 
_g/^"V 



(I _ g«-"V)-i 
g/^"V(I - g"^"V)^^ 



(18) 



It is clear from Eq. (18) that the computational effort in 
inverting M, and thus finding the inhomogeneous elec- 
tronic structure, is entirely determined by the complexity 
of inverting M"^". We also obtain the useful result that 
the electronic structure in the far field is easily deter- 
mined once the electronic structure in the near field is 
known. 

The local density of states in the near field is merely 



-^^ImG<,(x,x;a;) = -hm{2g{x,x;uj) + J dx' [[(M"^")-i(x,x';a;)] g(x', x; w)] 

-J dx' [[(M"-")(-i)*(x,x';-a;)] g*(x', x; -a;)] ^ J 



(19) 



while in the far field it is 



— — y^ImGo-(x, x;ti;) = ——lm\2g{x,x;w) 

a 

+ J dx'dx" [g(x,x';c^)V(x') [(M"-")-i(x',x";u;)]g(x",x;u;)]^^ 

-J dx'dx" [g*(x,x';-a;)V(x') [(M"^")(-i)*(x', x"; -a;)] g*(x", x; -a;)] ^ J. 



(20) 
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In the limit that the size of the near-field approaches a point, Eq. (20) is 

--y'lmGa(x,x;w) = --Im\2g{x,x;u>) + [g(x,0;a))Veff(w)g(0,x;w)]ii - [g*(x, 0; -w)V;fi-(-w)g*(0, x; -w)]22 }• 

a 

(21) 

The inhomogcneous part of the right-hand side of Eq. (21) has the following form if the off-diagonal elements of 
Veff(w) can be ignored (this is appropriate when self-consistency is not important): 



--Im {([Veff (a;)]ii - [Veff *(-a;)]22) (^(x, 0, u;)g{0, x, u) - /(x, 0, a;)/(0, x, to))} . 
n 



(22) 



Hence the spatially-dependent features reported in 
Ref. 26 are qualitatively retained in our more complete 
model. For a magnetic potential in the Born approxima- 
tion there would be no such signal, since 



Veff(a;) = 



Vs 
Vs 



(23) 



so the distortions in the LDOS of the two spin directions 
are equal and opposite. 

The local density of states under ideal conditions is di- 
rectly measured by a scanning tunneling microscope. The 
differential conductivity measured at a point x, voltage 
V and temperature T can be related to the local density 
of states at the tip location as follows: 



dV 



du)- 



9nsTM(w) /ImG'^(x,x; w) 



(24) 



Here e is the charge of the electron and nsTui'^j) is the 
Fermi function in the STM tip 



nsTM(w) 



1 + exp 



eV 



(25) 



The local density of states is proportional to the imagi- 
nary part of the retarded Green's function fully dressed 
by the interaction of the electronic system with the im- 
purity. 

Certain features of the above equations simplify their 
numerical implementation. Each angular momentum 
channel constitutes an independent block-diagonal sub- 
matrix in M(a;). Since the bare Green's functions in 
Eqs. (15) have analytic expansions in spherical harmon- 
ics, M(ti)) can be calculated quickly. These expansions 
are detailed in the Appendix. 

The value of R is governed by the longest-range poten- 
tial. In this paper that will be (5A(x). We model Vs and 
Vo with Gaussians of radius kp^ and evaluate Eq. (12) at 
T = 0. The (5A(x) for Vs = nNo\ J dxVs{x)\ = 1.75 and 
vq = ttNo I / d,xVo (x) I = is shown in Figure 1 . While 
oscillating with wavelength ~ nkp^, (5A(x) falls off to a 
negligible potential within lOA;^^. A typical radial grid 



I 

of 100 points provides a numerically robust solution. The 
self-consistent solution depends on the value of 



(26) 



which for niobium is 3.6 x 10~^ (Ref. 64). This is the 
single dimensionless parameter required to parametrize 
the frcK^-clcctron model of a siiperconductor. 

In contrast to the Kostcr-Slater technique described 
here, which exploits a physical distinction between the 
near field and the far field to accelerate the numerical 
calculation, the BdG equations treat the near-field and 
the far-field on the same level. A large R is desired for 
decent spectral resolution, but the possible size of R is 
limited by numerical constraints. Hence a numerical im- 
plementation of the BdG equations, by comparison, is 
typically slower and substantially less accurate than the 
Koster-Slater technique. 



< 

to 




r (units of kp 



FIG. 1. Change in the local order parameter, JA(x), for 
a magnetic potential strength Vg = No J drVs(r) = 1.75 as 
a function of the distance from the impurity r. The change 
becomes negligible beyond 10k . 



III. ANALYTIC SOLUTION OF THE POINT 
POTENTIAL 

A. Magnetic and Non-Magnetic Point Potentials 

Approximating the local potential by a delta function 

Vs + Vo 



V(x) = VS{x) = 







Vs~Vo 
leads to a simple expression for Eq. (4), 



6{x) (27) 
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l-g{0,0:u;)iVs + Vo) -/(0,0 
-r{0,0;iv)iVs + Vo) l + 5*(0,0; 



M"^"(a;)G(0,0;a;) = g(0,0;a;). 



(28) 



In principle M.{uj) can be found from the Green's func- 
tions in Eq. (15), however there is a divergence in the real 
part of ^(r; w) as r ^ 0. This divergence is coped with in 
Ref. 8 by discarding the divergent piece. This approxi- 
mation is essentially an assumption of strict particle-hole 
symmetry (not merely linearizing e(fc) around ep)- We 
now discuss the effects of this approximation on the lo- 
cal properties of the system in the norm,al state. The 
heuristics will be simpler for the normal state, but the 
conclusions also apply to the superconducting state. 



B. Particle-Hole Symmetry in the Normal State 

In order to focus on the spin-dependent potential, the 
non-magnetic potential will be set to zero. The nor- 
mal state properties can be obtained from Eq. (15) for 
Ao = 0. That yields a Green's function appropriate for 
an outgoing wave: 



Sf(x,x';u;) 



ttN, 



2_ ikpr 

e , 



(29) 



where uj is considered close to the Fermi surface so the 
change in momentum due to 7^ is negligible. The 
inhomogeneous local density of states for a delta-function 
potential using this Green's function is unphysical, since 
it diverges at r = 0. Now we examine the local density 
of states when the divergent real part of g is ignored. 



--ImGa{yi,x;to) = No I 1 



l-inNoVs)^ 



sin kpT 
(30) 



This expression yields the pathological behavior that the 
local density of states near a spin-dependent potential is 
exactly the same for spin-up electrons as for spin-down 
electrons. A spatial response in the LDOS to a spin- 
dependent potential that is identical for up and down spin 
electrons only occurs for certain band structures with 
at special energies {e.g. half- filled tight-binding models). 
The local density of states for Vs = 'itNo\Vs\ = 0.1 is 
plotted in Figure 2 as a function of distance from the 
potential. Although it is somewhat distorted from its 
homogeneous value, it does not show the spin- depen- 
dent asymmetry of a more realistic potential. An exact 
calculation for a Gaussian of range kp^ is also shown in 
Figure 2 for comparison. 

A more realistic approach to coping with the diver- 
gence in Eq. (29) without yielding the pathological re- 
sult of Eq. (30) is to average the real part of g over a 
range given by the assumed range of the potential. That 
yields a finite value for the local density of states at the 
potential, but does not control the behavior for small r. 
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To perform that task, we consider a "muffin-tin" Green's 
function. This function has the form 



5r(x,x';w) 



■kN, 



kpr 
= —ittNo — Tra 



r > Ro 

r<Ro (31) 



where a is the average of the real part of g over the range 
of the potential, and Ro is chosen so that the spatially- 
integrated spectral weight of the Green's function is un- 
changed (as required by probability conservation). We 
show in Figure 2 the local density of states calculated 
with this Green's function. In particular, the asymmetry 
in the response of spin-up electrons and spin-down elec- 
trons to a spin-dependent delta-function potential is gov- 
erned by this phenomenological parameter a. The agree- 
ment with the exact solution is good at the origin and 
far from the impurity. The muffin-tin Green's function is 
discontinuous, unfortunately, but yields a better approx- 
imation of the response of the system than the particle- 
hole symmetric approximation (discarding the real part 
of the Green's function). 

It is also possible to generate an asymmetry between 
spin-up and spin-down electrons by adding a non- 



1.15 I ---i. 

^pin up 



0.85 




Gaussian, range=kp' 
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FIG. 2. Spectral weight in the normal state around a mag- 
netic impurity with Vs = 0.1 as a function of the distance from 
the impurity r. The solid lines are exact solutions for spin-up 
(attracted by the impurity) and spin-down (repelled by the 
impurity) electrons for a Gaussian potential of range kp^ in 
a metal with a free-olcctron dispersion relation. The dashed 
line, which is the same for both spin-up and spin-down elec- 
trons, is the calculated spectral weight for the particle-hole 
symmetric delta-function potential model. The dot-dashed 
line is the result for a delta-function potential calculated us- 
ing the muffin-tin Green's functions. The muffin-tin Green's 
functions fix the pathological result of the symmetric model 
that the spectral weight is the same for spin-up and spin-down 
electrons. The muffin-tin parameter, a = 0.704, is determined 
by the range of the Gaussian potential and is therefore not a 
free parameter. 
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magnetic potential vq with the Vs to parametrize the im- 
purity, but still maintaining a particle-hole symmetric 
band structure. Since the response of a particle-hole sym- 
metric system to a potential does not depend on the sign 
of that potential, the vo is required to distinguish between 
electrons and holes. An electron feels a potential Vg + vq, 
whereas a hole feels a potential Vs — vq. This additional 
non-magnetic potential only breaks particle-hole symme- 
try locally (within the range of the potential), whereas 
for a realistic band structure the particle-hole symmetry 
is broken everywhere in the solid. 

For a Gaussian potential with range a, 



E 



1 



-(2n-3)!! 



(32) 



For the Gaussian potentials numerically calculated in this 
paper, a = kp^, so a = 0.704. 



C. Self-consistency within the Analytic Model 

As seen from Figure 1, the distortion of the order pa- 
rameter is short-ranged around an impurity. We may 
then consider the effect of the order-parameter distortion 
on the electronic structure to be parametrized by an effec- 
tive delta- function potential 5A(5(x) similarly motivated 
to the delta-function potentials for the magnetic and non- 
magnetic potentials. The potential V(x) is changed in 
the following way: 



V(x) = V5(x) 



Vs + Vo 5A \ 
6A Vs-Vor^''^- 



The relative effect of the 5A compared to the other 
two potentials is likely to be small for the potentials 
considered in this paper. Typically NoVs/kp ^ 1 or 
NoVo/kp ~ 1, and for niobium NgAo/kp = 3.6 x 
10~^. Even for a small coherence length of ^ = lOA;^^, 
NoAo/kp = 1.6 X 10~^. For convenience we define 
Sq = wNoSA. 



<?(0,0;u;) 



-ttNo a 



(35) 



The anomalous Green's function is given by Eq. (15), 
since it does not have a divergence problem as r — > 0, 



f{0,0;io) = -TrNo 



A, 



(36) 



This analytic model only has localized states in the £ = 
angular momentum channel, as expected for a delta- 
function potential. Those energies are 



where 



vJo± [(2^;A)^-4(i;g+7^)(^o-7^)] 



1/2 



Ao 

(37) 



7 = [(1 + a''){vl - vl - 5l) - 2avo - l] /2. (38) 

Eq. (37) reduces to a result obtained by Shiba^ when 

Vq = a = So = 0, a result obtained by Rusinov'' when 
a = So = 0, a result obtained by Salkola, Balatsky and 
Schrieffer^^ when a = 0, and a result obtained by us^^ 
when Vq = So = 0. 

One solution of Eq. (37) is a spin-up quasiparticle and 
the other is a spin- down quasiparticle. There may be 
only one real solution to Eq. (37); then only one i = 
localized state exists around the impurity. This occurs 
for large Vg. When Vs ~ the localized states are due 
to order parameter suppression, and the energies of the 
two spin states are degenerate. This follows from time- 
reversal symmetry in the absence of a magnetic potential. 
For small Vg the two energies are split by an amount 



2vsSoAo 



(39) 



D. Energies of localized states in the Superconductor 



E. Spectral Weight Asymmetry in the Analytic 
Model 



The energies of the localized states of angular momen- 
tum i correspond to the positive energies co£ = \rt\, where 

detM"^"(fi) = 0, (34) 

and the solution is traced to the ^-channel block of M 
(see Appendix). For the analytic model, M"^"(a>) is 
the matrix shown in Eq. (28), where 17(0,0; w) in the 
superconducting state is constructed similarly to that of 
the normal state, 



A spin-up quasiparticle consists of amplitudes for a 
spin-up electron (electron in a spin- up state), and a spin- 
up hole (electron missing from a spin-down state) . There- 
fore the spectral weight of a spin-up localized state will 
be divided between an electron-like pole in the spin-up 
band at lo = lUq (with weight A] {r;LUo)) and a hole-like 
pole in the spin-down band at w = —cOo (with weight 
A|(r; — Wo))- These two types of excitation are inde- 
pendently resolvable by a scanning tunneling microscope 
since at positive sample voltage relative to the tip, the 
STM places electrons in the sample, whereas at negative 
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sample voltage the STM places holes in the sample. We define the energy of the pole in the spin-up band to be uj^ and 
in the spin-down band to be uji. Even though cUq is always positive, can be positive or negative, and uj^ = —(^i- 

The spatial structure of the spectral weights of the spin-up band and spin-down band components of the localized 
state are given by 



A^{r; Lj) = — S{iJ - uJa)- 



2vs 



(^0 - Vs)Al + {vq + v^)!^"^ + (^0 + Vs)a^{Al - lo^) 
A2 



+ 



2{V0 + Vs)a OJ^f^^ - (1 + a2)(^2 _ ^2) (^^(^2 _ ^2) + 



Al 



— Wjj) X 



r <Ro 



^ A2-a;2 

+ {Vs +Vq)\1 + {Va - Vq) 



VQ + {vl -vl){a 



= - a > cos^ (kpr) 



sin^ {kpr) 



+2{vs + vo) \ vo - Vs + [1 - a{vo - Vg 



WW 



sin(fcFr) COs(fcirT") 
r>Ro 



(40) 



The above expression is set up for use as a muffin-tin 
Green's function. The construction of such a Green's 
function requires that the spectral weight integrates to 
one. However, if Ro 0, the above expression can be 
integrated over all space, and still yields one-half for the 
spin-up band and one-half for the spin-down band to or- 
der Ao/ep. Hence the quasiparticle is one-half electron 
and one-half hole. This provides additional confidence in 
the above expressions. 

The frequency-integrated weight at the defect of the 
two types of excitation can be calculated within the an- 
alytic model, and is (for 6o = 0) 

^T(r = 0) ^ 1 + 2a{vo - v^) + {1 + a'^){vo - v^ f 
A^{v = 0) l + 2aivo + Vs) + {l + a^){vo+Vs)^ ^ ' 

For a spin-down quasiparticle there is an electron-like 
pole in the spin-down band and a hole-like pole in the 

spin-up band and the relative weight is still given by the 
expression above. This expression for vq = was re- 
ported in Ref. 37. 

In the following section these results are compared with 
the numerical calculations of properties in the supercon- 
ducting state and in the normal state. 



IV. COMPARISON WITH NUMERICAL 
RESULTS AND DISCUSSION 

A. Non-Magnetic Impurity 

Even strong non-magnetic impurities at moderate con- 
centrations will not suppress the critical temperature of 
a superconductor^. Nevertheless, it was recognized early 
on^^ that the local order parameter may be affected. In 
Ref. 25 the effect of the non- magnetic impurity was calcu- 
lated in the far field by modeling the impurity potential 
with a phase shift. The phase shifts were evaluated for 
two models: a spherical square- well potential and a delta- 
shell potential. Self-consistency was ignored by only fo- 
cusing on regions far from the impurity where the change 
in the order parameter is small compared to its homoge- 
neous value. The order parameter change due to the im- 
purity was found to oscillate with the Fermi wavelength, 
and decay as r~^. 

Figure 3 shows the spectral weight A{r; uj) at several 
frequencies above the energy gap near a strong non- 
magnetic impurity with a Gaussian potential of range 
kp^ and strength vo = 7/8, calculated self-consistently 
for a superconductor with S^kp = 449 (niobium). The 
spectral weights are suppressed to approximately 30% of 
their homogeneous value at the center of the potential. 
Only continuum states are shown since no localized states 
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1 Each curve is displaced vertically by 0.1 for visibility 1 
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FIG. 3. Comparison of the spectral weight in the super- 
conducting state at several frequencies to the spectral weight 
of the normal state and the spatial dependence of the or- 
der parameter. All of these curves are indistinguishable, but 
they have been displaced for better visibility. The spectral 
weights are normalized to their value at large distances from 
the impurity r, and the order parameter is normalized to the 
homogeneous order parameter Aq. This figure indicates that 
normal-state properties drive the spatially-dependent features 
of the superconductor's spectrum. 

were found for this potential. The curves showing the 
spectral weight have been displaced from each other so 
that they may be distinguished. Also shown displaced in 
Fig. 3 is the spectral weight in the normal state, normal- 
ized to the spectral weight in the homogeneous metal. 
All the quantities plotted in Fig. 3 are identical to the 
accuracy of the calculation. Fig. 3 is an illustration of 
a relationship between the spectral weight in the nor- 
mal state and the spectral weight in the superconducting 
state, 

^(r;a;) = ^^cM (42) 

where A„(r) is the spectral weight in the inhomoge- 
neous normal state for energies near the Fermi surface 
and Asc{io) is the homogeneous superconductor's spec- 
tral weight as a function of frequency. 2No is the normal 
state's spectral weight far from the impurity. This ex- 
pression is valid for small r and small ui, the regime of 
interest for STM on a superconductor. For w of order 
Ao, Eq. (42) is valid for r < ^. 

We further illustrate the relationship of Eq. (42) in Fig- 
ure 4(a), which shows the LDOS for this non- magnetic 
impurity as a function of voltage and position calculated 
from Eq. (24) with T — ^Ag/kB- This temperature 
corresponds to 2K for niobium. There is no change in 
the energy gap due to this non-magnetic potential. Fig- 
ure 4(a) shows that it is merely the local amplitude of the 
spectral weight which is reduced this would manifest 
itself in a locally reduced oscillator strength for an opti- 
cal transition, or the reduction in the tunneling current 
for an STM, which is directly proportional to the LDOS 



shown in Figure 4(a). Figure 4(b) is an identical calcu- 
lation for a shorter coherence length, ^ = Wkp^. There 
appears to be little difference, although here a localized 
state exists near the continuum. 

Figures 3 and 4(ab) show for the non-magnetic impu- 
rity that the normal-state electronic structure determines 
the spatial dependence of the superconductor's A{r;Lo) 
for all frequencies including near the energy gap. The 
potential strength of the impurity is orders of magnitude 
greater than Ag, and thus locally mixes in states far 




FIG. 4. Differential conductance (LDOS) as a function of 
voltage and position calculated around a non-magnetic im- 
purity with vo = 0.875. The spectrum is suppressed sub- 
stantially in the vicinity of the impurity. The temperature is 
Ao/7.5A:b, which for niobium corresponds to about 2K. (a) A 
coherence length appropriate for niobium, ^kp = 449. (b) A 
much shorter coherence length, ^kp = 10. 
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from the Fermi surface in the homogeneous metal. These 
states are required to construct probabihty densities 
which are suppressed by 70% near the impurity. The spa- 
tial structure of the spectral weight in the normal state 
is essentially identical to that seen in Fig. 3 over an en- 
ergy range around the Fermi surface which is orders of 
magnitude greater than Aq. 

Once the normal-state band structure has been dis- 
torted by the presence of the non-magnetic impurity, 
superconductivity is a small perturbation within a nar- 
rower frequency range. The formation of the gap in the 
single-particle excitation spectrum in the superconduct- 
ing state is characterized by the "mixing" of electron 
and hole amplitudes to form quasiparticles near the gap 
edge. These quasiparticles, therefore, are constructed 
from single-particle eigenstates of the metal which have 
already been strongly distorted by the impurity poten- 
tial. 

Equation (42) has important implications for spec- 
troscopy on a superconductor, for one of the procedures 
for normalizing spectra taken at different lateral positions 
on a superconducting surface is to assume that the LDOS 
at a particular voltage much larger than Aq is the same. 
This is an attempt to correct for possible changes in the 
tip-surface distance upon moving the tip laterally. A 
small change in the tip-surface distance can have a strong 
effect on the tunneling current. Unfortunately this pro- 
cedure will prevent an experiment from seeing changes in 
the LDOS due to a non-magnetic impurity, including the 
conductance oscillations described in Ref. 26. 

We now discuss the properties of the order parame- 
ter. A(r) is self-consistently determined, and is shown in 
Fig. 3 for small r to be identical in spatial structure to 
the normal-state spectral weight, 

A(r) ^ AJrl ^ A{r;u;) 

Since a non-magnetic potential repulsive to electrons at- 
tracts holes, and A(r) depends equally on electron and 
hole amplitudes, one might expect a non-magnetic po- 
tential to have little effect on the spatial dependence of 
the order parameter. However the allowable maximum 
spectral density of holes depends on the spectral density 
of the electron band where the holes reside, so if most 
electrons are excluded from the site, holes will be effec- 
tively excluded as well. To emphasize this point we note 
that the scaled anomalous spectral weight ImF(r, r; lu) is 
identical to the scaled ^(r; u)) for the frequencies shown 
in Fig. 3 (and for all relevant frequencies for the self- 
consistency equation Eq. (12)). Since Imi^(r,r;a;) is 
proportional to the product of electron and hole ampli- 
tudes, and A{r; uj) is proportional to the electron ampli- 
tude squared, the spatial structure of the electron and 
hole spectral weights must be similar. They are, since 
the normal-state spectral weight An{r;u)) is roughly fre- 
quency independent around the Fermi energy over an en- 
ergy range much greater than Aq. 



We now comment on the lack of localized states near 
the non-magnetic impurity for (^kp = 449 and the small 
binding energy for the quasiparticle for ^kp = 10. The 
suppression of the order parameter near the impurity 
may be considered to form an attractive off-diagonal po- 
tential which may bind quasiparticles. Localized states 
created by order-parameter suppression would be dou- 
bly degenerate, due to the two possible spin states (see 
Eq.(37)). Since the quasiparticle is half hole and half 
electron, if the electron part is attracted and the hole 
part is repelled, one might expect the effects of such a 
non-magnetic potential on the quasiparticle to cancel. 
However the binding energy of the localized state is an 
order of magnitude smaller {ujo = (1 — 2 x lO^**) Ao) than 
that found in Section IV. D for a suppressed order pa- 
rameter via pairing suppression. This may be explained 
by the well-known repulsive effect (quantum reflection) 
of a strong attractive potential on a quantum-mechanical 
particle. We find also that in the case of the magnetic 
impurity that a large enough non-magnetic potential of 
either sign will suppress the binding of a quasiparticle to 
the impurity. We note here that the ratio of the non- 
magnetic potential to the off-diagonal potential {SA{r)) 
is much larger {vo/NgAo ~ lO"') than the ratios of the 
non-magnetic potentials to the magnetic potentials con- 
sidered below. 



B. Magnetic Impurity 

Recently, we have presented calculations of the LDOS 
(and thus the differential conductivity in an STM exper- 
iment) in the vicinity of a magnetic impurity^^. These 
calculations indicated that the spatial structure of the 
electron amplitude of the localized state differed strongly 
from the hole amplitude of the localized state. A further 
result was that the spectrum should recover to the ho- 
mogeneous spectrum within a few atomic spacings. Sim- 
ilarly motivated calculations of the LDOS due to the 
£ = localized state have been presented since then"*^. 
although these calculations did not address the contin- 
uum LDOS. The two models used in Ref. 38 were (1) 
a (5-function model solved using particle-hole symmet- 
ric Green's functions, but not self-consistently, and (2) 
self-consistent calculations for a two-dimensional tight- 
binding s-wave superconductor within the BdG equa- 
tions. The first method can only model the normal-state 
properties properly for a particle-hole symmetric band 
structure, such as at the Van Hove singularity in a two- 
dimensional tight-binding band structure. The second 
method must contend with numerical finite-size effects, 
which make it difficult to calculate the continuum states. 
A result obtained from the first method which is only true 
for special band structures is that the spatial structure of 
the electron and hole components of the quasiparticle are 
the same. The authors of Ref. 38 did raise the possibility 
of an additional non-magnetic potential as a source 
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FIG. 5. Energies of the quasiparticle poles as a function 
of the magnetic potential strength Vs for the first two an- 
gular momentum states around the impurity. The spin-up 
and spin-down labels refer to the band that the excitation 
resides in — an excitation with negative energy (hole-like) 
in a spin-up band is a spin-down hole. The solid line cor- 
responds to (J^F = 10, the long-dashed line to ^kp ~ 449, 
the short dashed line to the symmetric model of Shiba, and 
the dot-dashed line to the result calculated with mufHn-tin 
Green's functions and a — 0.704. At a critical value of 
~ v*Q, the up poles cross to negative energies and the 
down poles cross to positive energies, indicating a change 
in the character of the ground state. The kink evident in 
the solid and long-dashed lines is real, and due to the dis- 
continuous (at T = 0) change in A(x) at t;*o. The unim- 
portance of self-consistency can be gauged by the small dif- 
ference between the short-coherence length result and the 
long-coherence length result. 



of electron-hole amplitude asymmetry in the spatial 

structure of the localized state. We found^^ and will 
explore below that there is, for realistic band structures, 
electron-hole asymmetry without a non-magnetic poten- 
tial. 

We will begin with a discussion of the energies of the 
localized states aroiind a magnetic impurity and the spin 
character of those states. Solutions to Eq. (34) can be 
evaluated numerically. Figure 5 shows the dependence 
of the energies of the; localized state poles for the first 
two angular momentum channels on the strength of the 
magnetic potential. Results for a short coherence length 
(^ — 10k p^) are shown (solid line) as well as results for 
a long coherence length (^ = 449A;^^). The localized 
quasiparticle state for small Vg is the spin state attracted 
to the classical spin, which we will label up (t)^^. As 
the potential strength increases, the excitation energy of 
each angular momentum state decreases. At some critical 
value w*£ the localized state becomes a spin-down exci- 
tation, the energy changes abruptly, and then increases 
with increasing Vg. This behavior can be extracted from 
the analytic model (Eq. (37)) as well. 

Also shown in Figure 5 are the analytic results for the 
pole energies for a = (Ref. 8) and a = 0.704. The 



muffin-tin model is no better than the particle-hole sym- 
metric model in predicting the localized state energies. 
The muffin-tin model will prove more successful at pre- 
dicting the asymmetry between the electron and hole am- 
plitudes of the localized state. 

The unimportance of self-consistency for determining 
the energies of the localized states can be gauged by the 
small difference between the short-coherence length re- 
sult and the long-coherence length result. The most im- 
portant feature it determines in Fig. 5 is the size of the 
discontinuity in the localized state energy at v^q. This 
discontinuity is due to a discontinuous change in the or- 
der parameter at this magnetic potential strength, a re- 
sult pointed out in Ref. 37 which will be discussed more 
below. 

A magnetic impurity affects each member of a pair 
oppositely. An up-spin electron is attracted to the im- 
purity and a down-spin electron is repelled. The differ- 
ence in energy between the two time-reversed electron 
states leads to a locally suppressed pairing. Therefore, 
the energy needed to break a pair and create a localized 
quasiparticle with angular momentum £ in the vicinity of 
the magnetic impurity, when Vs < Vsi*, is reduced from 
2Ao to Ao -f LU£. One member of the broken pair is a de- 
localized spin-down quasiparticle at the gap edge (with 
energy Aq). The other member of the broken pair is a lo- 
calized spin- up quasiparticle with energy uj^. For Vg > v*^ 
there is a spin-up quasiparticle with angular momentum 
£ in the ground state of the superconductor^^. An es- 
sential point about the ground state of a superconductor 
containing classical magnetic impurities is that at T = 
when Vg < Vgi* (for all £) the ground state is composed 
entirely of paired electrons. When Vg > Vg(* for any £ 
then the T = ground state contains localized quasi- 
particles as well as pairs and there are new low-energy 
excitations including the re-formation of a pair as well 
as the excitation of a localized quasiparticle to a higher 
energy localized or continuum state. 

Evident from Eq. (5) is that the electron pole involves 
entirely single-particle states within the spin-up band for 
Vg < v*£, but involves entirely states within the spin- 
down band for Vg > v*^. The hole-like pole always in- 
volves single-particle states in the opposite spin band 
from the electron-like pole. The source of the quasiparti- 
cle amplitude for the various poles is indicated in Fig. 5. 

Figures 6(a), 6(b), and 6(c) show LDOS results for 
Vg = 0.52, 0.875, and 1.75 respectively. They show the 
state split off from the continuum, with a larger electron- 
like amplitude than hole-like amplitude (6(a)), and then 
lower in energy with an increased electron/hole asym- 
metry (6(b)). Finally the larger peak becomes hole- like 
(6(c)). The asymmetry between the electron-like and 
hole-like peaks becomes more pronounced as Vg increases. 
We note that the larger peak is always associated with 
the spin-up band, whereas the other is associated with 
the spin-down band. Despite the apparent differences 
in peak size, the spatially integrated electron spectral 
weight of the quasiparticle is equal to the spatially inte- 
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FIG. 6. Differential conductance (LDOS) calculated around a magnetic impurity with (a) Va = 0.5, ^/sf = 449, (b) Vs — 0.875, 
^kp = 449, (c) Vs = 1.75, ^kp = 449, (d) Vs = 0.875, ^kp = 10. All are calculated with fcsT = Ao/7.5. In the progression from 
(a) to (c) the asymmetry between the two peaks increases and the higher peak moves to lower energies, eventually becoming 
hole- like. 



grated hole spectral weight. The localized quasiparticle 

is always half electron and half hole for all potentials 
examined here. For Vs < v*q the spin-up band amplitude 
is electron-like and the spin-down band amplitude is hole- 
like. At v*Q (1.32 for niobium), the spin-up component 
becomes hole-like and the spin-down component becomes 
el(!(;tron-like, as required by the change in the spin of the 
excitation. 

Figure 6(d) shows the LDOS for a markedly different 
coherence length, ^ = lOkp^, and Vs = 0.875. It is evalu- 
ated for the same values of Ao/ksT as Figs. 6(abc) and 
looks almost identical to Fig. 6(b). Since the localized 
state is broadened by temperature through Eq. (24), this 



is a manifestation of the proportionality of the spectral 
weight to NqAo (Eq. (40)). Figure 7 shows the spectral 
weight at the origin for the £ = state and for the £ = 1 
state at its first maximum for Vs = 0.875 as a function of 
the inverse of the coherence length, which is proportional 
to NqAo- Figure 8 shows the spectral weight for ^ = 
Mdkp^ as a function of Vs for the spin- up and spin-down 
poles at the origin for the £ = state and at the first 
maximum for the £ = 1 state. It is clear that a non- 
magnetic potential is not necessary to obtain an electron- 
hole asymmetry. 

In Figure 9 we show the asymmetry at the impurity as 
a function of Vg for two values of ^ — a long coherence 
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FIG. 7. Spectral weight at the impurity (r = 0) for the first 
angular momentum channel, i = 0, and at the first maximum 
for the second angular momentum channel, £ = 1, for poles 
in both the spin-up and spin-down bands, as a function of 
the inverse coherence length, showing a linear behavior. The 
magnetic potential strength is Vs = 0.875. 

length appropriate for niobium, and a short coherence 
length. From Figure 7 it should be evident that the 
asymmetry is not sensitive to ^. It is, however, pre- 
dicted extremely well by the normal-state spin-up and 
spin-down band spectral weight asymmetry at the impu- 
rity (also shown in Figure 9). We can therefore conclude 
that as with the non-magnetic impurity (Figure 4), the 
spatial structure of the spectral weight is a normal-state 
property. We further show in Figure 10 for Vg = 0.875 
the £ = 0, spin-up band and spin-down band projections 
of the normal-state spectral weight to compare with the 
localized state spin-up band and spin-down band spec- 
tral weights for two values of ^fc^- The normal-state and 
long- coherence length calculation are practically indis- 



For the short- coherence length calculation the 



o 

03 



3 



'CD 
g 
CO 



0.0002 

0.0001 

0.0000 
0.00004 



■G 0.00002 



o 

03 



0.00000 



■ Spin up 

Spin down 






/=0 










r Spin up 






; Spin down 







0.0 0.5 1.0 1.5 

Magnetic Potential Strengtli (v^) 

FIG. 8. Spectral weight at the impurity (r = 0) for the 
£ = channel for ^kp = 449 as a function of magnetic po- 
tential strength for poles in both the spin-up and spin-down 
bands. The spectral weight in the spin-up band pole of the 
£ = localized state saturates at large Vs ■ Also shown are the 
spectral weights at the first maximum for the £ = 1 localized 
states. 
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FIG. 9. Ratio of the spectral weight in the spin-up band 
and in the spin-down band at the impurity (r = 0) as a func- 
tion of magnetic potential strength. This is plotted for the 
normal state £ = projected Green's functions (solid line) 
as well as for the localized states for niobium {^kp = 449, 
long dashed line), for ^kr = 10 (dot-dashed line), and for the 
muffin-tin model. The muffin-tin model is only successful for 
Vs < 0.5, but that is duo to a breakdown in describing the 
normal state. The normal-state electronic structure is a good 
predictor of the superconductor's electronic structure for the 
entire range of Vs- 

effect of an exponential envelope is also visible. In the 
analytic result th e exponential envelope should have a 
range R = 7r^/2\/l — (wo/Aq)^, which corrects to better 
than 1% the discrepancy in Figure 10. The power-law 
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FIG. 10. Spectral weights for £ — localized state in the up 
and down bands for Va — 0.875. The solid line is the normal 
state £ = projected spectral weight, the long-dashed line 
is the localized state in a superconductor with ^kp = 449, 
and the dot-dashed line is for ^kp = 10. The inset shows 
the spectral weight multiplied by r'^ to remove the algebraic 
decay. The normal-state and long-coherence length results 
are practically indistinguishable. The deviation shown in the 
short-coherence length superconductor's spectral weight is fit 
to within 1% by the exponential decay factor described in the 
text. Hence the spatial structure of the spectral weight of 
the superconductor's localized state is well-predicted by the 
normal state spectral weight. 
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fall-off and exponential envelope can be seen directly 
from Eq. (15) and Eq. (40). 

We can summarize these comments with a general 
equation, similar in concept to that for the non-magnetic 
impurity, Eq. (42). That is, for a localized quasiparticle 
state of spin a' , the spectral weight of a localized state 
with angular momentum t would be 



(44) 

where i? is a normalization factor so that the spectral 
weight of the state integrates to one, and Ano-{ic,£) is 
the angular-momentum £ projection of the normaZ-state 
spectral weight in the spin a band. We note that for 
small r there is an approximate relationship between the 
superconducting state's spectral weight and the normal 
state's spectral weight in each spin band. 



1 

— J dLjA„{r;u) = Ana{r), 



(45) 



where Aq <C -E <^ ep- This, in connection with Eq. (44), 

implies a dependence on the normal-state structure of the 
continuum spectral weight around the magnetic impurity. 

We now return to the structure of A(x). This quan- 
tity, which is not directly observable, formed the focus of 
several investigations of the local structure around a 
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FIG. 11. Order parameters as a function of distance from 
the impurity r calculated for several magnetic potential 
strengths for (a) S^kp = 449 and (b) S^kp = 10. In both 
cases there is a discontinuous change in the order parameter 
when Vs passes through the critical strength v*q. 
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FIG. 12. Order parameters at the impurity (r = 0) as a 
function of magnetic potential strength Vs for two values of 

the coherence length. The discontinuity in the order param- 
eter at v*Q is much larger for the short-coherence-length su- 
perconductor. 

magnetic impurity. The oscillation of the order parame- 
ter around a magnetic impurity was first evaluated with- 
out self-consistency^^~^^. A self-consistent calculation of 
the order parameter at the impurity and very far away 
for weak impurity potentials was done by Schlottmann^^. 

As shown in Figure 1, for large values oi Vg, A(x = 
0) < 0. Sign changes in A, as seen in pair tunneling, 
have been suggested for magnetic impurities in the bar- 
riers of Josephson junctions^^"^*. The sign change in 
A(0) occurs (at T = 0) precisely at v*q. Due to the spin 
and frequency symmetries of Eqs. (4-11), the anomalous 
spectral weight Imi^(r, r, cj) associated with the spin-up 
pole is always equal and opposite to the anomalous spec- 
tral weight associated with the spin-down pole. As the 
pole in the spin-up band goes from electron-like {u > 0) 
to hole-like (lu < 0) and the pole in the spin-down band 
goes from hole-like to electron-like the contribution to 
A(0) changes sign abruptly at T = 0. The A(r) resulting 
from several values of Vg and two values of the coherence 
length are shown in Figure ll(ab). The discontinuity 
at v*Q is more pronounced for shorter coherence lengths 
since the localized state's spectral weight is more con- 
centrated at the impurity (Eq. (44)). A(0) as a function 
of Vs is shown in Figure 12 for two values of the coher- 
ence length. At T > the transition would be smoothed 
somewhat. 

The behavior of A(0) as a function of Vg comes from 
the introduction at w*o of a quasiparticle into the ground 
state of the system. The spin-up quasiparticle localized 
near the impurity in the ground state suppresses the lo- 
cal order parameter. For time-reversal invariant poten- 
tials one cannot make A(r) negative by inserting a single 
quasiparticle, since the suppression from one quasiparti- 
cle is cancelled by the lack of suppression from its un- 
excited Kramers doublet partner. For a spin-dependent 
potential, however, the anomalous spectral weight near 
the impurity may be almost entirely contributed by the 
single low-energy localized state. When a quasiparticle is 
present in the ground state, the ground state has spin 5 
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^65.38 



up""'"" and a negative A(r)^'. Exciting the low-energy 
state for Vg > ij*q removes the spin-up quasiparticle, and 
therefore increases A(0), whereas excitation of quasi- 
particles typically reduces A(x) (which is the case for 
Vg < v*q). Also, exciting the low-energy state reduces 
the induced spin of the superconductor at the impurity. 



C. Combined Magnetic and Non-Magnetic 
Potentials 

We now discuss the addition of a non-magnetic poten- 
tial to the magnetic potential. It had been suggested^^ 
that introducing a vo with a Vs would provide electron- 
hole asymmetry. We find that it does change the asym- 
metry, which we show in Figure 13 for a particular Vs, but 
that once again this is a normal-state property. The rela- 
tionship between the normal-state spectral weights and 
the superconducting-state spectral weights of Eq. (44) 
still holds. Introducing vq also alters the localized-state 
energies (see Eq. (37)), which we show in Figure 14 for 
Vs = 0.875 and ^kp = 449. The presence of a non- 
magnetic potential may affect the value of flo'"*^- We 
show in Figure 15 a partial diagram of the ground state as 
a function of the parameters Vs and vq ioi ^kp = 449. We 
note that the boundary between the two ground states is 
not shifted much from (,kp ~ 449 to ^fc^? = 10, hence the 
condensation energy is not very significant in determin- 
ing this boundary. This fact, and the observation that 
the energy of the localized state (in the absence of self- 
consistency) relative to the homogeneous order parame- 
ter does not depend on any parameter of the supercon- 
ducting state (Eq. (37)) suggests that an understanding 
of the transition at v*q may be assisted by normal-state 
concepts. We hope to address this issue in a future 
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FIG. 13. Ratio of the spectral weight at the impurity in 
the spin-up band to the spin-down band for Vs = 0.875 and 
^kp = 449 as the non-magnetic potential vo varies. In a 
similar result to that seen in Fig. 9, the normal-state spectral 
weights are good predictors of the superconducting state's 
spectral weight. We note tiiat the curves are not symmetrical 
around vo = 0, which results from a realistic band structure 
without particle-hole symmetry. 
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FIG. 14. The energy of the spin-up pole is shown as a func- 
tion of non-magnetic potential strength for Vs = 0.875 and 

^kp = 449. The energy of the spin-down pole is just the 
negative of the energy of the spin-up pole. The analytic mod- 
els do not perform particularly well in reproducing the pole 
energies, although the muffin-tin model does show a similar 
asymmetry around vo = to the numerical calculations. 

publication. 

We now mention an example designed to showcase 
the drawbacks of Ginzburg-Landau theory. We present 
in Figure 16 the order parameter for a mixed magnetic 
and non-magnetic impurity {Vg = 7/8, vq = —7/4, and 
^kp = 449). The order parameter is everywhere larger 
than in the homogeneous superconductor, however, the 
presence of a localized state within the gap indicates that 
superconductivity has been weakened around the impu- 
rity. Since Ginzburg-Landau theory focuses on the or- 
der parameter and ignores the quasiparticle structure, a 
Ginzburg-Landau perspective would incorrectly predict 
an enhancement of superconductivity in the region. 

We conclude by attempting to provide some guidance 
for attempts to extract impurity potentials from STM 
measurements. Figure 17 shows normal-state dl/dV's for 
various potentials. These curves should also represent 
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FIG. 15. Calculated boundary between two ground states 
around the magnetic impurity for ^kp = 449. For a large 
enough magnetic impurity strength a quasiparticle is bound in 
the ground state. The minimum magnetic impurity strength 
depends on the non-magnetic impurity strength. For still 
larger magnetic impurity strengths there would be ground 
states with more than one bound quasiparticle. 
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FIG. 16. The change in order parameter around an attrac- 
tive non-magnetic potential combined with a magnetic po- 
tential for ^kp = 449. The order parameter is larger at the 
impurity than in the homogeneous superconductor, yet there 
exists a localized state in the gap. 

the frequency-averaged spectral weight measured in the 
superconducting state (see Eq. (45))^®. 

The enhancement or suppression of spectral weight 
near the origin is particularly sensitive to vq. A measure- 
ment of this quantity, the energy of the localized state 
and the asymmetry of the electron and hole amplitudes at 
the impurity ovcrconstrains vq and Ug, given an assump- 
tion of the shape of the potential. To extract informa- 
tion about the potential's detailed shape would require 
a fitting procedure using the differential conductivity at 
various positions. If, for some reason, the spin-down am- 
plitude were too small to measure, it may remain possible 
to constrain the potential strength using the frequency- 
integrated spectral weight and the localizcd-state energy. 

It seems appropriate to mention again the tendency 
to normalize spectra according to the LDOS measured 
at energies much larger than Aq. Since thc^ normal- 
state LDOS near the impurity changes substantially in 
the presence of magnetic or non-magnetic impurities, an 
experiment performed using such a normalization pro- 
cedure would yield impurity parameters of questionable 
validity. 
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FIG. 17. Spectral weights in the normal state as functions 
of the distance from the impurity for several combinations 
of magnetic and non-magnetic potentials (a) only a magnetic 
potential, and (b) a fixed magnetic potential Vs = 0.875 and 
a varying non-magnetic potential. By making measurements 
around the impurity in the normal state (or integrating the 
superconductor's spectrum over a frequency much larger than 
Ao), information about the structure of the impurity may be 
obtained. 



D. Pairing Suppression 

The pairing potential, 7(x) in Eq. (12), may also have 
spatial structure. When this parameter is changed it in- 
duces a change in the order parameter which produces 
an off-diagonal potential felt by the quasiparticles. We 
set Vs = Vq = so that there is no magnetic or non- 
magnetic potential to compete with the order parameter 
change. Figure 18 shows the order parameter around a 
short-range suppression, 



7(x) = [l - e-^*'^'''' 



Jo, 



(46) 



for two values of the coherence length. The order param- 
eter is strongly suppressed and since 7(0) = 0, A(0) = 0. 
For long coherence lengths, this change in the order pa- 
rameter has no effect on the local density of states (shown 
in Figure 19(a)). It is possible to localize quasiparticle 
states, however, at shorter coherence lengths. These 



can produce features in the local density of states which 
are visible. One such case is shown in Figure 19(b). The 
energy of the localized state is Wo = (1 — 4 x 10~^)Ao. 
Whereas a non-magnetic potential changes the local 




r (units of l<p ) 

FIG. 18. Order parameter as a function of distance r from a 
defect with a suppressed pair potential, but no single-particle 
potential. 
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FIG. 19. Differential conductance (LDOS) around a defect with a suppressed pair potential, but no single-particle potential, 
(a) ^kp = 449, fesT — Ao/7.5 — there is no evidence of any change in the spectrum duo to the order parameter suppression 
shown in Fig. 18. (b) ^kp = 10, ksT = Ao/lOO — a localized state very near the continuum enhances the continuum edge 
seen in tunneling near the defect. 

density of states without significantly changing the energy gap, a pairing suppression has a very weak effect on both, 
especially in the long- coherence length limit. 

V. STRONG-COUPLING AND ANISOTROPIC ORDER PARAMETERS 

A few observations are in order concerning the extension of this formalism to systems where the homogeneous 
order parameter has important frequency or spatial structure. The Gor'kov equation (Eq. (4)) changes due to the 
more general form for the off-diagonal potential originating from the order parameter. Taking this opportunity to 
generalize, 

J dx"dx"' [5(x - x")(5(x" - x'") - g(x, x"; a;)V(x", x'"; uj)] G(x"', x'; w) = g(x, x'; uj) (47) 

where 

\[x,x,u>)-\^ ,5A(x",x"';a;) Vhiix";io)6{x" - x'") J ' ^^^> 



The diagonal terms are general potentials, possibly 
frequency-dependent, effective on spin-up electrons (K-f) 
and spin- up holes {Vh-\) for u> > 0. Since this poten- 
tial is diagonal in frequency, as is the Gor'kov equation, 
the frequency structure of the order parameter does not 
add any additional complication to numerically solving 
the Gor'kov equation. However, the addition of spatial 
structure to the pairing has added another integral over 
the volume to the Gor'kov equation, and thus dramati- 
cally increased the size of the matrix which needs to be 
inverted. Fortunately, the range of the order parameter 
in |x — x'l is truncated by the range of the pairing inter- 
action. For the isotropic-gap superconductor we have 



I 

considered for most of this paper, the effective pairing 
interaction is modeled by a delta function in space. It 
is possible to obtain anisotropic order parameters, in- 
cluding the d-wave order parameters possibly appropriate 
for high-temperature superconductors, merely by consid- 
ering pairing with nearest neighbors on a square tight- 
binding lattice. This model has been implemented with 
the BdG formalism for non-magnetic impurities'^^"'^^ and 
magnetic impurities^^ and vortices™ in rf-wave supercon- 
ductors. When formulated on a lattice the addition of 
nearest- neighbor pairing multiplies the rank of the matrix 
M"^"(a;) hy 1 + z, where z is the coordination number 
of the lattice. 
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The order parameter's frequency dependence compli- 
cates the self-consistency equation (Eq. (12)). It must 
now be solved for each frequency: 

A(r,r';w)Z(r,r';u;) 

/OO 
den{e) [ReF{r, r'; e)] {K+{r, r'; e, to) - U,) 
-OO 

u[l-Z{r,r';iv)] 

/OO 
den{e) [ReG(r,r';e)] Js:_(r,r';e,w), (49) 
-co 

where K± are kernals of the pairing interaction and are 
different for each mechanism of superconductivity. They 
can be determined from the homogeneous solution. Uc 
is a Coulomb factor. Z, the quasiparticle weight, is 
solved for sclf-consistently. Incorporating these strong 
coupling effects allows a determination of the effect of 
the frequency-dependence of the pairing interaction on 
the electronic structure around a defect. While selecting 
the particular location where the STM tip is similar to 
selecting the momentum of the quasiparticles of interest, 
selecting the STM voltage indicates which order param- 
eter frequency one wishes to probe. 



the inhomogcneous normal metal multiplied by a decay- 
ing exponential whose range is determined by the energy 
of the localized state. Since the spin-up band LDOS in 
the normal state differs from the spin-down band LDOS 
in the normal state, the electron-like pole of the localized 
quasiparticle will have different spatial structure than the 
hole-like pole of the localized quasiparticle. 

The self-consistent calculations described here have 
been performed with a new, powerful, Koster-Slatcr tech- 
nique which allows the Gor'kov equation to be solved in 
principle exactly. Although we have only presented cal- 
culations for weak-coupling isotropic order parameters 
within a free-electron model, we have formulated the ex- 
tension of the new technique to strong-coupling pairing 
and general band structures and order parameter sym- 
metries. 
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VI. SUMMARY 



APPENDIX 



The local electronic structure around an defect reflects 
both the properties of the defect and the medium it is 
embedded in. For a strong non-magnetic or magnetic im- 
purity in a superconductor, the distortion of the normal- 
state properties by the strong impurity plays a vital role 
in the response of the superconducting medium. The 
LDOS for the inhomogcneous superconductor can be re- 
lated to the LDOS for the inhomogcneous normal metal 
via equations hke Eq. (42) and Eq. (44). In the case of 
a non-magnetic impurity with no localized states around 
it, the LDOS for the inhomogcneous superconductor is 
merely the normalized LDOS for the inhomogeneous nor- 
mal metal multiplied by the homogeneous superconduc- 
tor's density of states. This should suggest some caution 
for the method used to normalize STM spectra taken at 
different places on a superconductor's surface. 

For the case of a localized state (such as around a mag- 
netic impurity) with angular momentum £, the LDOS for 
the state is the angular-momentum-^-projected LDOS of 
I 



The expansion of the Green's hmctions of Eq. (15) 
suitable for a three-dimensional spherically symmetric 
situation are detailed here. The homogeneous Green's 
functions depend on r and r' through r, r', and cos 7 = 
(r • r')/rr'. Then the Green's functions can be expanded 
in Legendre polynominals P^(cos7). 



5f(r,r';a;)= 

g{r, r',cos7;a;)= 



2i+l 

4-77 



J29e{r,r';u;)Pe{cos^) (50) 



and 



gi{r,r'\uj)=2TT I Pe{x)g{r,r' ,x;lo). 



(51) 



Evaluating Eq. (51) for both the normal and anomalous 
Green's functions of Eq. (15) yields 



ge{r, r';u) 



fe{r,r';uj) = - 



TT 



UJ 



irr' — 1 



J,+ . ({1 - ^f:s^l|i}r<)llfl^{{\ - ^^;^/ar>) 



(52) 



%i ({1 + ^/j^\li}r<)E'lU{^ + ^/J^\li}r>) 
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(53) 



where J^, H^^\ and Hf are standard Bessel functions, 
(r^) is the smaller (larger) of r and r', w is in units 
of Ao and r and r' are in units of kp^. The Green's 
functions arc in units of Ng. 

The Gor'kov equation, Eq. (4), can now be written in 
a form diagonal in 

Gi{r,r'-u:)= (54) 

g£(r,r';w)+ / r^dr„g£(r,r„; w)V(r„)G£(r„,r';w). 
Jo 



r(2) 



Thus the three-dimensional integral has been reduced to 
a one-dimensional radial integral. Since the numerical in- 
version procedure depends on inverting a matrix whose 
rank is proportional to the number of spatial points con- 
sidered, this reduction to a one-dimensional integral dra- 
matically increases the speed of this calculation over a 
calculation for a three-dimensional potential which is not 
spherically symmetric. 
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